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Abstract 

Observing prices of European put and call options, we calibrate exponential Levy models 
nonparametrically. We discuss the implementation of the spectral estimation procedures for 
Levy models of finite jump activity as well as for self-decomposable Levy models and improve 
these methods. Confidence intervals are constructed for the estimators in the finite activity 
case. They allow inference on the behavior of the parameters when the option prices are 
observed in a sequence of trading days. We compare the performance of the procedures for 
finite and infinite jump activity based on real option data. 
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1 Introduction 



Exponential Levy models are frequently used for the purpose of pricing and hedging. They take 
jumps of the price process into account and they allow to model heavy tails in the returns appro- 
priately. Moreover, they are capable of reproducing not only the volatility smile but also the fact 
that it becomes more pronounced for shorter maturities. While recovering these stylized facts, the 
structure of exponential Levy models is at the same time easy enough to allow robust calibration 
procedures. 



The calibration has mainly focused on parametric models, cf. Barndorff-Nielsen (1998); Eber 



lein, Keller, and Prause (1998); Carr, Geman, Madan, and Yor (2002) and the references therein. 



First nonparametric calibration procedures for finite activity Levy models were proposed by |Cont 
and Tankov ( 2004b ) as well as Belomestny and Reif5 ( 2006a ) . In these approaches no parametriza- 
tion is assumed and thus the model misspecification is reduced. Recently, the spectral calibration 
procedure of the latter authors was studied further in two directions. The asymptotic confidence 



sets constructed by Sohl (2012) allow a more accurate statistical inference in the exponential Levy 



model while the method of Trabs (2011) extends the spectral calibration to the infinite activity 



case, more precisely to self-decomposable Levy processes. 

The aim of this work is to improve the spectral calibration method and to apply it to real 
option data. In contrast to the method of Cont and Tankov ( 2004b ) the spectral calibration is a 
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straight forward algorithm, where no minimization problem has to be solved. Therefore, already 
the originally proposed methods by Belomestny and Reifi (2006b) and by Trabs (2011) are quit 
fast owing to the FFT-algorithm. Whereas previous works focus on the asymptotic theory, we 
concentrate on the application of the method to realistic sample sizes. By a preprocessing of the 
data and a different choice of the weight functions we reduce the loss of the estimators significantly. 
These considerations are demonstrated by several simulations from the model by ^Merton ( 1976 1 



and from the variance gamma model, introduced by Madan and Seneta ( 1990 ) and Madan, Carr. 



and Chang (1998). 



Exponential Levy models are studied in a wide range of pricing problems, for instance by 



Asmussen, Avram, and Pistorius (2004); Cont and Voltchkova (2005); Ivanov (2007) and the 



references therein. For their application, it is necessary to control the statistical errors in the cali- 
bration procedures, see Cont (2006). In general, there are two types of errors. The misspecification 
is the deviation from the model, while the calibration error is the deviation within the model. The 
approach of nonparametric calibration has the advantage that it reduces the error due to mis- 
specification. The remaining calibration error can then be assessed by means of the confidence 
intervals. Nonparametric confidence intervals and bands for Levy densities have been constructed 
by iFigueroa-Lopez (20111 based on high frequency observations. Sohl (2012) derived asymptotic 
confidence sets for the calibration of the risk neutral measure and based on observations of option 
prices and not on historical data. However, they turn out to be too conservative in simulations. 
To describe the behavior of the estimators more precisely, we construct confidence intervals based 
on a finite sample variance. As demonstrated by simulations these intervals perform well in terms 
of size and coverage probabilities. 

We use real data of vanilla options on the German DAX index to compare the finite activity 
model to the self-decomposable one. Both models achieve a good calibration in the sense that the 
residuals between the given data and the calibrated model are small. In view of studies of |Carr| 



et al. ( 2002 ) and Ait-Sahalia and Jacod ( 2009 1 , which indicate that a pure jump model should have 



a Blumenthal-Getoor index which is positive or even greater than one, this is surprising since the 
finite variation self-decomposable model has a Blumenthal-Getoor index equal to zero. Applying 
the calibration to a sequence of trading days, we obtain the evolution of the model parameters in 
time. 

This paper is organized as follows: In Section 2 we state the model and the general estima- 
tion method. This is made precise for the finite activity case and the self-decomposable case in 
Sections 3 and 4, respectively. Simulations are contained therein. In Section 5 the confidence 
intervals are constructed and their performance is assessed in simulations. In Section 6 we apply 
the methods to real data and discuss our results. We conclude in Section[7] The appendix contains 
the proof of a corollary and a detailed documentation of our implementation. 



2 Model and estimation principle 

Denoting the riskless interest rate by r > and the initial value by 5*0 > 0, the price of a stock is 
given by 

St = 5'oe''*+^* 

where Xt is a Levy process with characteristic triplet (ct^,7, such that e'^* is a martingale. 
Furthermore, we assume the jump part of Xt to be of finite variation. Therefore, £[6^^*] = l,t > 0, 
implies the martingale condition 

Y+7 + /_"(e"-l)Kdx)=0. (1) 

So far, nonparametric calibration methods exist in two different setups: 



(FA) Xt has a finite activity, that is A := < oOj with an absolutely continuous jump measure 

dCont and Tankovl |2004b| IBelomestny and Reifil |2006al ISohll |2012|. 



(SD) Xt is self-decomposable with a = 0. Especially, ly can be characterized by a k-function 
via i/(da;) = fc(a;)/|a;| dx,x € M \ {0}, where k is increasing on ]R_ and decreasing on M_|_. 
Additionally, a := k{0+) + fc(O-) < oo is assumed ( |Trabs[ [20Tl| ). 

Throughout we measure the time t in years. Typical parametric submodels of (FA) and (SD) are 
given by Examples [T] and [2] respectively. We will use them to study the performance of estimation 
methods in simulations. 



Example 1 (Merton model). Merton ( 1976 ) introduced the first exponential Levy model. Therein, 
the jumps are normally distributed with intensity A > 0: 



u{x) 



X 



■ exp 



X e 



A realistic choice of the parameters is rj ~ —0.1, v — 0.2 and A = 5. Together with the volatility 
cr = 0.1 this determines the drift 7 = 0.379 using the martingale condition Q. 

Example 2 (Variance gamma model). Let [Wt] be a standard Brownian motion and (Gt) an 
independent Gamma process with mean rate one and variance rate v that is Gt ^ T{t/i>, l/v). 
Madan and Seneta ( 1990| ) defined the variance gamma process with parameters a, v and Q as the 



time changed Brownian motion with drift Xt — OGt + aWa^,t > 0. The characteristic function 
and the k-Function of (Xt) are given by 

(pt{u) = (1 + i6vu + G^vi? jiy^^" and 



(x) 



4a;<0} 



{x) + -e l{^>o}(a;), u,x G 



with rip/m '■= ■\/ 0-2/4 + 0-2 ± 9v/2, respectively. In our simulations we use the parameters 
cr = 1.2, V = 0.2 and 9 = —0.15. The value of 7 = 0.141 is again given by the martingale condition. 
These choices imply a — kvG[^+) + fcyG(0— ) = 10. 

Let us fix a maturity T > 0, define the negative log~moneyness x := \og{K / So)~rT and denote 



call and put prices by C{x,T) = 5'oE[(e'^^ 
In terms of the option function 



e=")+] and V(x,T) = 5oE[(e"= - e^^)+\, respectively. 



0{x) :-- 



'S^'C{X,T), 
So'V{x,T), 



x>0, 
x < 0, 



our observations are given by 



Oj =Oixj) + S,ej, j = l, 



(2) 



with noise levels Sj > and independent, centered errors ej, satisfying Var(ej) = 1 as well as 
supjE[ej] < 00. The observation errors are due to the bid-ask spread and other market frictions. 
To estimate the Levy triplet, we apply the Levy- Khintchine represen tation of the characteristic 
function ipriu) ■= E[e™^^]. The pricing formula of Carr and Madan (1999^ 



and the martingale condition ([T]) yield 
1 



e™^0(a;) dx 



- (Pt{u- i) 
u{u — i) 



T 



\og{l + iu{l + iu)FO{u)) 



+ i(cr^ +7)m+ / (e 



(3) 



l)e''v(x)dx. (4) 



Interpolating {xj, Oj)j=i,...,Af, we obtain empirical versions O and ip of the option function and the 



characteristic exponent, respectively. While the theoretical results (Belomestny and Reifi 2006a 



3 



Trabs 2011) concentrate on a linear interpolation of the observation, an additional smoothing by 



using B-splines of degree two might improve the estimators. In Section [3^2] we provide simulations 
with both interpolation methods to investigate the practical influence. 

Given ip, we can estimate the characteristics of the process from the spectral representation. 
Regularization of the procedure is achieved by cutting off frequencies larger than the regularization 
parameter U > 0. Since (FA) and (SD) need to be considered separately, the precise estimators 
are given in the Sections [3] and |4] Note that in both cases correction steps are necessary to satisfy 
the shape restriction of v and the martingale condition ([ij (see Appendix |B . 2| for details). If the 
latter one would be violated, the right-hand side of the pricing formula ([3| could have a singularity 
at zero and thus we could not apply the inverse Fourier transform to obtain an option function 
from the calibration. 

A critical question is the choice of the regularization parameter U. To circumvent the problem 
in simulations, we use an oracle cut-off value, that is the risk minimizing U. To calibrate real 
data, we employ a simple least squares approach. From theoretical consideration a penalty term. 



as used by Belomestny and Reifi (2006b), is necessary to avoid an over-fitting. Nevertheless, our 



practical experience with this method shows that the above mentioned correction steps, which 
are not included in the theory, lead to an auto-penalization: Owing to the Fourier techniques, 
a rigorous fit leads to high fluctuations of the estimator of the nonparametric part and thus the 
correction has a significant effect which in turn worsens the fit. Therefore, the least squares choice 
of tuning parameter works well at least for small noise levels. Additionally, an upper bound for 
the cut-off value excludes estimations with a too high variance (cf. Section [s]). 



3 The finite activity case 
3.1 The estimators 

In the (FA) framework we deduce from Q the identity 



-u^ + i{a'^ + -i)u + ((7^2 + 7 - A) + J"/x(u) with /i(a;) e^z^(x). 



The estimators of the parameters are defined by Belomestny and Reifi (2006a) as follows: 



:= 



Ke{i]:{u))w^ {u)Au 



u 



7 := —(P 



Im('(/;(u))w^(M)dM, 



Re{'4)(u))w^{u)Au, 



(5) 
(6) 
(7) 



■u 



with suitable weight functions w^, w}^ and w^. We propose to choose the weight functions 
differently to the weight functions used by Belomestny and Reifi (2006b). The idea is that the 



noise is particularly high in the high frequencies and thus it is desirable to assign less weight to the 
high frequencies. A smooth transition of the weight functions to zero at the cut off-value improves 
the numerical results significantly. Therefore, we would like the weight function and its first two 
derivatives to be zero at the cut-off value. With the side conditions on the weight functions this 
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leads to the following polynomials: 

2s / „ s 2s+2 / „ s 2s+4 



-4(2. + 7)(-) +(2. + 9)(-) ), 

^"A^(-) I ((2^ + 3) (^) - 4(2^ + 5) (^) + 6(2. + 7) (^) ' 
_4(2,s + 9)(^)'^''%(2. + ll)(^)' 



where all three functions are extended continuously by zero outside [—U,U], c„,c^ and cx are 
normalization constants and s reflects the a priori knowledge about the smoothness of v. The gain 
of the new weight functions is discussed in Section [3?2| Using the shifted 



to estimate the jump density, we apply the idea from [Belomestny and Reifi (2006b) owing to a 



better numerical behavior. Therefore, we define the estimator 



{x) (8) 



with the empirical version tp,^ of -01/ and a flat top kernel whose support is [—U^U]: 

(1, |u|<0.05, 
cxp( --P(7ff7)r^)"^ ), 0.05<|u|<l, (9) 



0, \u\ > 1. 



3.2 Simulations I 



Let us first describe the setting of all of our simulations. In view of the higher concentration 
of European options at the money, the design points {xi, . . . ,xn} are sampled deterministically 
from a normal distribution with mean zero and variance 1/2. The observations Oj are computed 
from the characteristic function tpT using the fast Fourier transform. The additive noise consists 
of independent, normal and centered random variables with variance |rO(a;j)p for some relative 
noise level r > 0. By choosing the sample size N and the deviation parameter r, we determine the 
noise level of the observations. According to the existing theoretical results, it is well measured 
by the quantity 

e := A^/^ + A^/^||(5||;=o with A max {x.j - Xj„i), 

J —2, ... ,7V 

which takes the interpolation error and the stochastic error into account. The interest rate and 
time to maturity are set to r = 0.06 and T = 0.25, respectively. 

Using the Merton model with the parameters of Example [l] we investigate the practical influ- 
ence of two aspects of the procedure, which are mentioned above. The interpolation of the data 
{xj,Oj) with linear B-splines is compared to the use of quadratic B-splines. The latter prepro- 
cessing is an additional smoothing of the data which achieves significant gains for higher noise 
levels. The other point of interest is the choice of the weight functions. Since it is known from the 
theory that the noise affects mainly the high frequencies, the polynomial weight functions greatly 
reduce the variance of the estimator. These improvements are illustrated in Figure [T] In the case 
of a we calculate the RMSE from 500 Monte-Carlo iterations with and without quadratic splines 
and polynomial weight functions, respectively. This is done for different noise levels, whereby r 
decreases from 0.03 to 0.015 and N increases from 50 to 400, simultaneously. 
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Basic estimator 
Quadratic B-splines 
Polynomial weight 
Both 



Figure 1: RMSE of a for different noise levels with 500 Monte-Carlo iterations in each case. 
Usage of the linear and quadratic spline interpolation as well as usage of the weight function of 



Belomestny and Reifi ( 2006a| ) and the polynomial weight. 



4 The self— decomposable framework 
4.1 Construction of the estimators 

Recall that cr = is assumed in the (SD) setting. While the Blumenthal-Getoor index is zero 
in this case the parameter a describes the degree of activity of the process on a finer scale. The 
exponential scaled k-function h(.(x) := sgn(x)e^A:(a;), x S K, which is the counterpart of [i in the 
(FA) case, will be the natural object to estimate. To calibrate the self-decomposable model, we 
need a different representation of ^ than before because of the infinite activity of these processes. 
Trabs| (poTH Prop. 3.2) showed for u 7^ 



V'(u) = D(u) + i'^u — a log(|u|) + 



U3 



+ p{u) — iju -\- i / J^ke{v)dv, 



(10) 



where s is the smoothness of k away from zero, aj := fc(-'')(0+) + fc(^)(0-),j = 1, . . . ,s, the function 
D is constant on the real half lines and the remainder p satisfies ||m*^~^p(u)||oo < 00. Owing to the 
polynomial decay of p, estimators of 7 and a can be defined analogously to Section [3] 



lm{ip{u))wll (u) du, 



a 



-u 
u 



Reixjj{u))w^{u) du. 



The weight functions are chosen such that they filter the coefficients of interest (cf. [Trabs 2011 



Ass. 3) . Choosing and as polynomials, these integral conditions lead to a system of linear 
equations which determines the coefficients in the polynomials. Since the smoothness s is not 
known to the practitioner, the weights should satisfy the condition for some upper bound Smax- 
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In practice the choice 



^max — 6 and 



135135 ( 
1536 ' 



451*5 + 560u^ - 1890u^ + 2376u" - lOOlu"), 



<(u) = C/-X(^), 



1024 



4729725m*' - 56756700u« + 187297110u' 



231891660u^2 + 96621525u") 



works fine. The estimation of the nonparametric object fcg relies on the equation ip' = ij + HFke 



which fohows from (10). Hence, we need an empirical version of ij)' , too. Instead of the estimator 



proposed by Trabs 



20111, we define 



ke{x) := 



V-i [(_^ _ i^'{u))FWk{u/U)\ (a:), x > 0, 
J"-^ [(-7 - ii^'{u))FWk{-ulU)\ (a;), a: < 



with a one-sided kernel function Wk that satisfies 
Assumption 1. We assume 

(i) suppIVfc C [0,oo), 

(ii) / W^fe = 1, / x'Wkix) dx = /or / - 1, . 



- 1 and x'^'-^Wk{x) e 



(in) J \u\'^'^°'^'^\J-Wk{u)\ < oo for some upper bound a of the true a. 

Since we know the position of the jump of fcg, the application of a one-side kernel function 
allows to estimate the k-function on the whole real line. The asymptotic analysis of Trabs (2011 1 
carries over to this estimator as it is shown by the following corollary. Its proof is postponed to 
Appendix [A| Recall the definitions of the Sobolev-type class GsiR, a) and HsiR, a) from Def. 4.1 
and Def. 4.4 in |Trabs| ( [20TT| ). 

Corol lary 1. Assume s > l,R,a > and the conditions on A and 5 from Thm. 4-5 w | Trab^ 
'2011). Furthermore, let 7 be an estimator ofj satisfying supp E-p [I7— 7p]^^^ ^ g(2s+i)/(2s+2'j'a+5) _ 



Using the cut-off value U :— e 2/(2s+2Ta+5)^ obtain 



sup E-p[||/Ce 



^^||2^jl/2 < g2s/(2s+2Ta+5) 



Remark. For a < (5s — 1) / (2T) the above defined estimator 7 satisfies the assumed asymptotic risk 
bound of Thm. 4.2 in |Trabs| ( |2011| | and thus the corollary, restricted to "P G GsiR, ") ^1 TisiR, a), 
is applicable to the proposed procedure. 

The condition (iii) in Assumption [T] is a smoothness condition on Wk. To be adaptive in a, 
we construct the kernel as Wk{x) = P{x)w^{2x — 1) G C°°(K) with a polynomial P{x) of degree 
TO + 1. The coefficients of P are given by an (m + 1) dimensional system of linear equations, 
defined by Property (ii), which can be solved numerically. Rearrangement of fee ensures the 
necessary monotonicity of a k-function. 



4.2 Simulations II 

We simulate the variance gamma model from Example[2]and with the observation setting described 
in Section 3.2 The behavior of the proposed method for different noise regimes is shown in Table[l] 
As one expects from the theory, 7 can be estimated better than a. In the realistic setting of 100 
observations with one percent noise, the error of 7 and a lay below three and seven percent, 
respectively. The one-sided kernel achieves good results for kf.{x) even for small a; S M, indicated 
by Figure [2j However, there are still problems at zero. This emphasizes the importance of an 
separate estimator of a. 
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N 


50 


100 


100 




T 


0.01 


0.01 


0.05 


RMSE 


'7(RMSE) 


0.0014 


0.0004 


0.0074 


a (RMSE) 


0.4196 


0.6669 


2.0231 


RMISE 


(RMISE) 


0.5930 


0.4980 


0.4754 



Table 1: RMSE and RMISE using 1000 Monte-Carlo simulations of the variance gamma model. 




Figure 2: Estimated functions (left) and k (right) in a simulation of the variance gamma model 
with N = 100 and r = 0.01. 

5 Confidence intervals 



Sohl| ( 20T2| shows asymptotic normality of the estimators in the (FA) setup. These result may be 



used to construct confidence intervals. By (2.8) in Sohl (20121 it holds 



Aa^ - ^2 ^ f Re{J'fi{Uu))wl{u)du +^ f Re{Ai>{Uu))wl{u)du (11) 

'-^ Jo U Jq 

with Alp = ip — ip. The first term is the deterministic approximation error and the second term is 
the stochastic error. The choice of the cut-off value U allows a trade-off between these two errors. 
In order to construct confidence intervals, the cut-off value U is chosen such that the approximation 
error is asymptotically negligible. The term A'0([/u) is a logarithm, whose linearization at one 
we call Cn,u{u)- We denote by TZn.u{u) the remainder of this linearization. The stochastic error 
may then be decomposed as 

Re{Ai,{Uu))wl{u)du ^ ^ Re{CN,u{u))wl{u)du + ^ Re{TZN,u{u))wl{u)du. 

With an appropriate choice of the cut-off value the second part of the stochastic error is asymp- 
totically negligible. We will call the first part in the decomposition the linearized stochastic error. 
Confidence intervals may be constructed in two different ways. One can derive confidence intervals 
either from the asymptotic variance or from the finite sample variance of the linearized stochastic 
errors. We will follow the second approach. Nevertheless, the confidence intervals are asymptotic 
in the sense that the remainder term of the stochastic errors and the approximation error are 
considered as negligible. 

We assume that the noise levels of the observations ([2| are given by the values 6j — 5{xj), 
j — 1,...,N, of some function 5 : K — M+. The observation points are assumed to be the 
quantiles xj = H^^{j/ [n + 1)), j = 1, . . . , A^, of a distribution with a c.d.f. : M — > [0, 1] and 
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p.d.f. h. For the definition of the confidence intervals we need the generahzed noise level 

g{x) = S{x)/y/h(^), (12) 

which incorporates the noise of the observations as well as the density of the observations. Instead 
of assuming that the observation points are given by the quantiles of h one may also assume that 
the observation points are sampled randomly from the density h. 

For we calculate the finite sample variance s^2 of the linearized stochastic errors using 
z^)/2 and (6.17) from [S^ ( |2012| 

21 




CN,uiu)wl{u)du 



1 ^ 2- 

CN,u{u)wl{u)du 



47r 



Re 



Jo 



fuiu)fu{v)T-^ig{y + eY){U{u - «))e^^'^'("'+''')/2dudz; 



fu{u)fu{v)F-\g{y + ef){U{u - ^;))e^-'^'("'+"')/2d«d^;) , (13) 



+ Re 

where 9 — T{a^ + 7) and fu{u) — wl.(u){—u'^ + iu/U) exp(— rj^/i([/u)). The confidence interval 
for A is based on the similar finite sample variance of the corresponding linearized stochastic error. 
The only difference in the calculation for 7 is that we use Im(z)^ = Re(|zp — •z^)/2. We observe 
that i){x) in Q involves ipi, instead of ^. Thus the confidence intervals for vlx) are based on the 
linearization t/(u) of Aipi, = -ipv — 'i'l^- The variance s^(j.) of the linearized stochastic errors is 




Re(/:5V,£/(u)e~"^")wi(u)du 



^0 



f gu{u)gu{v)J''\e'^yg{y + T-/f){U{u^v))e^^'^'^^'+'''^/^dudv^ , (14) 



with gu{u) := w],{u){v? +iu/U) cxp{—TJ-i'{Uu)~ixUu). The confidence interval for ^{O) is based 
on both the linearization £^ jj multiplied by the weight function wjj, as well as the linearization 

£jv^(7 multiplied by J^^ v'^wj^{v)dv/2 — w\ J^^ w^(v)dw, cf. the definition of Wq by S 
We denote by Sa-2, s^, s\, Sj,(j.) and s^(o) the estimated standard deviations which are o 
by substituting cr, 7, A, and v by their respective estimators in (13), (14) and the analogous 
expressions. This yields feasible confidence intervals. For a level a > the (l-Q;)-confidence 
intervals are then given by 



(2012). 
btamed 



la- 
If 
h 



= [7 - S^qa/2,1 + S^qa/2] 

= [A - Sxqa/2, A + Sxqa/2], 

= [£>(a;) - Sr,{x)qa/2,i'{x) + Si.(x)9a/2], 

= [£>(a;) - Sry{0)qa/2,l'{x) + S^(0)9a/2], 



(15) 



where Qa denotes the (1 — a)-quantile of the standard normal distribution. 
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Figure 3: True Levy density (blue, solid)^ point- 
wise 95% confidence intervals (red, dotted) and 
100 estimated Levy densities (grey) from a 
Monte Carlo simulation. 



Figure 4: Boxplots of i>(0) and i>(0.4) 
from 1000 Monte Carlo iterations with 
true value (blue), theoretical quantiles of 
the linearized stochastic error for 25% and 
75% (red, solid) as well as for 2.5% and 
97.5% (red, dashed). 



We examine the performance of the confidence intervals by simulations from the Merton model 
with parameters as in Example [ij with interest rate r = 0.06 and with maturity T — 0.25 as in 
Section [3^ We sample = 100 strike prices deterministically and take the relative noise level to 
be r = 0.01. To coincide with the theory, we interpolate the corresponding European call prices 
linearly and use the weight functions by Belomestny and Reifi ( 2006b ) with smoothness parameter 
s — 2. In the real data application in Section 6 take advantage of the above improvements. Upon 
fixing an upper bound of 28 for the cut-ofF values, we perform 1000 Monte Carlo iterations. 

Figure [3] illustrates the true Levy density, pointwise 95% confidence intervals and the first one 
hundred estimated Levy densities from the Monte Carlo simulation. The graph shows that the 
confidence intervals describe well the deviation of the estimated jump densities. The negative 
bias around zero might come from the smoothing which naturally tends to smooth out peaks, cf. 

Chap. 5.3). The confidence interval Iu{a) is larger than indicated in Figure [3j By 



(Hardle 1990 



the smoothness of the estimated curves the larger variance of !)(0) leads to an increased variance 
of i'{x) in a neighborhood of zero. This is the second reason why deviation of the estimated curves 
is larger at zero. Figure [l] shows boxplots of i>(0) and 1^(0.4), the true value, the theoretical 25% 
and 75% quantiles as well as the theoretical 2.5% and 97.2% quantiles based on the linearized 
stochastic errors. We see that j)(0) has a negative bias and that the distribution of i){OA) fits well 
to the theoretical quantiles of the linearized stochastic errors. 



We asses the performance of the confidence intervals (15) with levels a = 0.5 and a = 0.05 
for the parameters a^, 7, A, z/(0) and ^{OA) in a Monte Carlo simulation with 1000 iterations. 
We approximate the coverage probabilities of the confidence sets by the percentage of confidence 
intervals which contain the true value. Table [2] gives the approximate coverage probabilities from 
the Monte Carlo simulation. 
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a' 


7 


A 




v{QA) 


a = 0.5 


54% 


47% 


50% 


49% 


47% 


a = 0.05 


97% 


95% 


91% 


98% 


95% 



Table 2: Approximate coverage probabilities of (1 — a)-confidence intervals from a Monte Carlo 
simulation with 1000 iterations. 
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Figure 5: Number of observed prices of put and call options during May ,2008. 



6 Empirical study 

We apply the calibration methods to a data set from the Deutsche Borse database EurexQ It 
consists of settlement prices of European put and call options on the DAX index from May 2008. 
Therefore, the prices are observed before the latest financial crises and thus the market activity is 
stable. The interest rate r is chosen for each maturity separately according to the put-call parity 
at the respective strike prices. The time to maturity T is measured in years. The number of our 
observations is given in Figure [5] 



6.1 Comparison of (FA) and (SD) 

Let us first focus on one (arbitrarily chosen) day. Hence, we calibrate the option prices of May 
29, 2008, with all four different maturities to both, the (FA) and the (SD) setting. The results are 
summarized in Table|3]and Figure[6j Using the complete estimation of the models, we generate the 
corresponding option functions O. They are graphically compared to the given data points and we 
calculate the residual sum of squares RSS — J^ji^j ~ ^i^j))'^- For all maturities both methods 
yield good fits to the data. However, for longer maturities, especially the calibration of options 
with four months to maturity, minor problems occur in the (SD) calibration. The calibration at 
other trading days confirms this weakness of the (SD) method for larger T. This coincides with 
the asymptotic analysis of Trabs (20111 and Corollary [I] where longer maturities lead to slower 
convergence rates of the risk. 

Moreover, Figure [g] shows that the estimated option function O which results from the (SD) 
calibration does not exactly recover the peak of C In all maturities and in both models the Levy 
density has more weight on the negative half line and thus there are more negative jumps than 
positive ones priced into the options. This coincides with the empirical findings in the literature 



(see eg, Cont and Tankov 2004a). 

In Table [3j we see a tendency that higher values of a in (FA) correspond to higher a in the 
(SD) model. The scatter plot in Figure [t] displays the pairs (a, a) in the estimations of all trading 

^provided through the SFB 649 "Economic Risk" 
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N 


61 


55 


101 


106 




T 


0.136 


0.233 


0.311 


0.564 




a 


0.109 


0.111 


0.106 


0.097 


(FA) 


7 
A 


0.224 
3.498 


0.163 
1.716 


0.176 
1.862 


0.193 
2.314 




y/RSS 


0.003 


0.007 


0.005 


0.007 




7 


0.337 


0.297 


0.378 


0.227 


(SD) 


a 


25.313 


30.084 


36.123 


16.557 




Vrss 


0.006 


0.006 


0.021 


0.014 



Table 3: Results of the estimations from option prices from May 29, 2008. 



days with maturity in July or August. Note that there are no more than 60 observations for each 
point and thus the stochastic error is relatively big. Nevertheless, the positive correlation between 
a and a is confirmed weakly and it is in line with the expectation that fluctuations of the stock are 
modeled by the diffusion part in the former model and by the infinitely high activity of the jump 
part in the latter one. The other way around, this relation shows that a is indeed a meaningful 
measure of the jump activity in the (SD) setting. 



6.2 (FA) across trading days 

The aim of this section is twofold. By considering more than one day we investigate the stability 
of the (FA) estimation procedure. Moreover, calibrating the model across the trading days in 
May, 2008, shows the development of the model along the time line and with small changes in the 
maturities. To profit from the higher observation number, we apply the calibration procedure for 
the (FA) case to the options with maturity in September and December. 

To apply the confidence intervals (15) of Section [sj we need the noise function g from (12). 
By a rule of thumb we assume S to be 1% of the observed prices 0{xj) (cf. Cont and Tankov 
2004a p. 439), for j = 1, . . . , A^. The density h of the strikes is estimated by a triangular kernel 
estimator, where the bandwidth is chosen by Silverman's rule of thumb. 

The estimations of the parameters are displayed in Figure [Sj Furthermore, the 95% confidence 
intervals for the December options are shown. The estimated volatility a fluctuates between 0.1 and 
0.12. The confidence sets imply that there is no significant difference of the two maturities. Both 
7 and A decrease for higher durations: On the one hand the curves of December lay significantly 
below the ones of September, on the other hand the graphs have a slight positive trend with 
respect to the time axis, which means with smaller time to maturity. 

Figure [9] displays the estimated jump densities. All jump measures have a similar shape 



which is in line with real data calibration of Belomestny and Reifi ( 2006b ) . In contrast to Cont 
and Tankov] ( |2004b I the densities are unimodal or have only minor additional modes in the tails, 
which may be artefacts of the spectral calibration method. The tails of D do not differ significantly, 
while the different heights reflect the development of the jump activities A. Showing the point- 



wise confidence sets around v^x) on a fine grid. Figure 10 confirms this suggestions. There is an 
obvious trend to small negative jumps in all data sets which is in line with the stylized facts of 
option pricing models. 



7 Conclusions 

To reduce the model misspecification it is reasonable to use a nonparametric model for option 
pricing. However, the nonlinear inverse problem, which occurs by calibrating the model, is more 
difficult to solve than parametric calibration problems and needs non-standard algorithms. We 
could improve the existing spectral calibration procedures for the finite activity (FA) Levy model 
and the self-decomposable (SD) Levy model. Owing to the fast Fourier transform, the method 
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Figure 7: Pairs {a, a) estimated from options from each market day in May, 2008, with Maturity 
in July (red diamonds) and August (green circles) and regression Hne (black). 



is computationaUy fast and admits convincing results in simulations and real data applications. 
Determining the finite sample variances of the linearized estimators, we obtain confidence sets, 
which allow a precise analysis of the noise of the estimators. 

Our empirical investigations show that both models can be calibrated well to European option 
prices. However, (FA) is more suitable for longer maturities. Using the derived confidence sets, 
we can observe significant changes of the (FA) model over time. While the volatility has no 
systematic trend, the jump activities decrease for longer maturities and thus the Levy densities 
become flatter. 

To avoid misspecification of the model, we are convinced that the nonparametric approach 
should be pushed forward theoretically and in practice, in particular, in view the high number 
of available observations in highly liquid markets. Of further interest would be extensions of the 
method to models whose jump part is not of finite variation as well as the application to hedging 
and risk management problems. 



A Proof of Corollary [T] 



W.l.o.g. it suffices to estimate the risk on M+. We decompose the risk to 



<3 



e|lL2(R+) 



L2(R+), 



{k,*{UWkiU<)){x)~k,{x) dx + 3Ep[|7-7|2] / \UWk{Ux)\^ dx 



+ 3E 
=:D + G + S. 



T-'{{i^'iu)-i;'{u))TW,{^)){x) 



dx 



The deterministic error term D can be bounded exactly as in the proof of (Trabs 2011 Thm 4.5) 



with r = 1 and use of the properties of the kernel Wt- The assumption on 7 and the choice of U 
yield G < [/e(4s+2)/(2s+2Ta+5) ^ ^4s/{2s+2Ta+5) _ n remains to estimate the stochastic error term 
S. Applying Plancherel's equality, the bound of \^p' — ip'\ from (Trabs, 2011 p. 27) and the decay 
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Figure 8: At each market day in May, 2008, estimated cr^ (top), 7 (center) and A (bottom) from 
options with maturities in September (dashed) and December (solid) and confidence intervals 
(dotted) for the latter ones. 





Figure 9: Estimation of v for maturity in September (left) and December {right). 
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Figure 10; Empirical confidence sets (dotted) of ('{x) (solid) in tlie estimation from option prices 
of May 6, 2008, with matm'ity in December and estimated jump densities from May 2 to 12 with 
the same maturity (grey). 



of TWk, we obtain 

^ , ||2 



S <Er[\\{^'iu)~^'iu))TWk{ 



U' 



< 



E[\HO - +E[|^(x(a - 0){x)){uf])\FWk{'^)\^ Au 



< ^2jj2Ta+5 _ g4s/(2s+2Ta+5) 

B Implementation in R 

B.l Global parameters 

We give here a detailed description of our implementation of the calibration method for the fi- 
nite activity case. Some global parameters to adjust the code and some auxiliary functions are 
necessary. These are namely Zeta needed for simulations, CFHat and ZetaHat to reproduce obser- 
vations from the estimated model, the Fourier transform function FT and a continuous logarithm 



logc. The calibration itself is done in the function calibration documented in Section B.2 We 
use these functions either to calibrate simulations or to estimate real data in Section iB.Sl 
There are a few parameters to control the execution of the program: 

• model permits the values "Merton" , "Kou" and "real data" . In the first two cases observa- 
tions will be simulated with the corresponding model. Picking the last value, the user has to 
supply the necessary data in the workspace of R before running the script (vector of strikes 
k, vector of option prices op, interest rate r and maturity T). 

• In case of simulations the user has to specify noiseLevel, the sampleSize and the number of 
Monte Carlo iterations in monteCarlo. Furthermore, design decides whether the observed 
strikes are sampled "deterministic" or "random" . 
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The choice of the cut-off parameter U is adjusted by mode, fn Belomestny and Reifi (2006b I 
are three possibihties proposed which are (in modifications) also available in this implemen- 
tation. Possible values for mode are: 

— "oracle": U* and U* minimize the discrepancy between the estimators and the given 
values, separately for (cr^,7, A) and i^. 

— "flat": The cut-offs correspond to points where the estimators stabilize: 



U* — argmin 



u 



dU 



dU 



aU 



a> 0, 



L2 



In our simulations and real data estimations we used a = 10^^. 
"PLS": The common cut-off U* solves 

N 



inf 

u 



Y,\C(K,;Tu)-Y,\^ + a f K{x)\^ dx 
.i=i •'^ 



a > 0, 



where C{K; Tjj) is the price at strike K computed using Levy triplet Tu = {(tu, lUi i'Uy )• 
In our simulations and real data estimations we use a = 0. 

— "fix": The same cut-off values Uf ix and Uf ixNu are used to estimate the parameters 
and the jump density, respectively. 

• If the variable linear is true then the observation will be interpolated linearly. Otherwise 
quadratic B-splines are used. 



• The smoothness of the jump density is set in parameter s (in the context of Belomestny and 



Reifi (2006a) s implies a smoothness s = 2*s of v). 



We remark that the library cobs is necessary to use this program. It provides the spline 
interpolation. 



B.2 The function calibration 

Given a vector of log-strikes sk and corresponding call option prices snop this function performs 
the calibration of the model and returns the cut-off values U,Ui,, the model parameters (7,7, A 
and the estimated jump density on a grid x. If a simulation is performed then also the quadratic 
L^-error of v is given back. 

c alibr at ion<— function ( sk , snop){ 
# return 

list (U=U, UNu=UNu, sigmaHat=sigmaHat , gammaHat=gammaHat , lambdaHat= 
lambdaHat , x=x , nullat=nullat , nuError2=nuError2 ) 

} 

To approximate the O-function from the given discrete points (sk, snop), we switch from the 
log-scale to the normal one, because we can use the convexity of the function K i-> C{K, T) = 
e~'""^E[(S'T — K)~^] for the quadratic spline interpolation. Furthermore, we extrapolate the obser- 
vations by adding boundary points whose position is set by extrapolation. To get an option 
value for the new strike near zero we use the slope of C{K, T) in K at 0. In the case of quadratic 
spline interpolation the function cobs from package cobs is used to interpolate and evaluate the 
result on the logarithmic scale on a fine grid. We get the points (knew, opnew) and use the put-call 
parity to have a discrete version of the function fc i— > ©(/c — rT). 
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extrapolation<— 0.005 



extraK<-rep (0 , length (sK) +2) 
extraK [ 1 ]<— cxt r apolat ion *sK [ 1 ] 
extraK [2 :( length (sK) +1) ]<-sK 

extraK [ length (sK) +2]<— sK [ length (sK) ]*l/extrapolation 

extraOp<-rep(0 , length (sK) +2) 
extraOp [ 1 ]<— 1— exp(— r*T) *extraK [1] 
extraOp [2 : ( length (sK) +1) ]<— snop 
extraOp [length (sK) +2]<-0 

BI<-log( extraK [1] ) 

BE<-log ( extraK [ length ( extraK ) ] ) 

knew<-seq(BI ,BE, length = 2'12) 

if ( linear ) { 

OpofK <— approxfun ( extraK , extraOp , method—' linear ", rule=2, ties = 
mean) 

opnew <— Opo£K(exp(knew) ) 
} else { 

css<— cobs ( extraK , cxtraOp , constraint—' convex" , nknots^min( 100 , length 

(sK)-2) ,degree=2) 
opnew<— predict ( CSS ,exp(knew)) [ ,2] 

} 

opnew<— opnew^max(0. ,1 — exp(knew— r*T) ) 
opnew<— pmax(0 , opnew) 

We can now apply the Fourier transform to get first the estimate of f h- > lpt{v — i) = 1 + + 
iv)FO{v) and calculate then 'il){v) = l/T\og((p{y — i)) by profiting from symmetry. Now we have 
(v, psi) and define L as the index where v has value 0. 

z<-FT(knew-r*T, opnew, TRUE) #x=k-rT , z=FO(v) 
v<— z$u 

phi<— z$fy *1 i *v*(l + l i *v)+l 
psi<— rep (0 , length ( phi ) ) 

psi [( length (v)/2) : length (v) ]<-l/T* logc ( phi [( length (v) /2) : length (v) ] ) 
psi [1 : ( length (v) /2-l)]<-Re( psi [( length (v) -1) : ( length (v) /2+1) ] )-li*Im 
(psi [( length (v)-l):( length (v)/2+l)]) 

L<-length(v) /2 

In the following we calculate for numerous cut off values U the estimators (7^,7 and A and save 
each in sigma2HatCur, gammaHatCur and leimbdaHatCur, respectively. The cut-off values are 
multiples of the mesh of v and their number is given by cutOff Iterations. Since the calculation 
time in the "PLS"-mode is much longer than in the other modes but the picked cut -off is smaller 
in general we choose less iterations in "PLS" . If the mode is "fix" , only one iteration with the fixed 
cut-off suffices. Because of a large bias for small U we start with (v[2]-v[l] )*11 and therefore 
end at (v [2] -v [1] )* (10+cutOff Iterations) . Later we need bestError for the decision in the 
oracle-mode. With use of symmetry it is enough to calculate (v_cut, psi_cut) on the positive 
half-axis. 

if (mode=" oracle" | mode="flat") 



18 



cutOfflte rations <— 130 
else if (mode="PLS" ) 

cutOfflteration s<— 9 
else if (modc^" fix" ) 

cut OffIterations<— 1 

signia2HatCur<— rep (0 , cutOfflterations ) 
gammaHatCur<— rep (0 , cutOfflterations) 
lambdaHatCur<— rep (0 , cutOfflterations ) 
bestError< — 1 

for ( i in 1 : cutOfflterations ) { 
if (mode=" fix" ){ 
UCur<-Ufix 

gridU<-round(Ufix/(v[2] -v [1] ) ) 
} else { 

UCur <- (v[2] -V [1] ) *( i+10) 
gridlJ <- i+10 

} 

v_cut<-v [L : (LfgridU) ] 

psi _cut<— psi [L : ( Lf gridU ) ] 

} 

Within this loop we approximate the integrals in ([5]) - ([7|) with a composite trapezoidal rule and 
use the polynomial weight functions from above. 

w<— rep ( 1 , length ( psi _cut ) ) 
w [ length ( psi _cut ) ]<— . 5 

wFkt<- f unction (x) { (2 *s + l)*x ^ (2 *s )-4* (2*s +3)*x ~ (2* s +2)+6* (2* s +5)*x 

^(2*s+4)-4*(2*s+7)*x^(2*s+6) + (2*s+9)*x'(2*s+8)} 
weight<— wFkt ( v_cut /UCur) 

weight [ length (weight ) ] <— weight [ length (weight ) ] —2*sum(w* weight ) 
sigma2HatCur [ i ]<—sum(w* weight *Re( psi _cut ) ) 

sigma2HatCur [ i ]< — 2*sigma2HatCur [ i ] /sum(w*abs ( v_cut ) "{2}*weight) 

wFkt<- f unction (x) {x" (2*s + l)-3*x" (2*s+3)+3*x" (2*s+5)-x~ (2*s + 7)} 
weight<— wFkt ( v_cut /UCur) 

gammaHatCur [ i ]<—suni(w* weight *Im( psi _cut ) ) 

gammaHatCur [ i ]<— gammaHatCur [ i ] /suni(w*weight*v_cut )— sigma2HatCur [ i ] 

wFkt<- f unction (x) { (2 *s+3)*x ' (2 *s )-4* (2*s +5)*x ~ (2* s +2)+6* (2* s +7)*x 

" (2*s+4)-4*(2*s+9)*x^(2*s+6) + (2*s + ll)*x'(2*s+8)} 
weight<— wFkt ( v_cut /UCur) 

weight [ length (weight ) ] <— weight [ length (weight ) ] —2*sum(w* weight *v_ 

cut"2)/v[L+gridU]"2 
lambdaHatCur [ i ]<— suni(w*Re( psi _cut )* weight ) 

lambdaHatCur [ i ]< — lambdaHatCur [ i ] /sum(w* weight )+gammaHatCur [ i] + 
sigma2HatCur [ i ] /2 

The choice of U depends on the mode. In "oracle" in every iteration step we calculate the distance 
between the estimators and the true values as an error term and select the minimizing U. In the 
case of "fix" the cut-off is unique. 

i f ( mode^" o r a c 1 e " ) { 
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err or<— l*abs ( sigma2HatCur [ i ] — sigma " 2)+l*abs (gammaHatCur [ i ] —gamma 

) +l*abs ( lambdaHatCur [ i ] — lambda ) 
if((i==l) I error <bestError ) { 
U<-UCur 

sigmaHat<— sqrt (pmax( , sigma2HatCur [ i ] ) ) 
gammaHat<— gammaHatCur [ i ] 
lambdaHat<— lambdaHatCur [ i ] 
bestErro r<— error 

} 

} 

if (mode=" fix" ){ 
IK-UCur 

sigmaHat<— sqrt (pmax(0 , sigma2HatCur [ i ] ) ) 
gammaHat<— gammaHatCur [ i ] 
lambdaHat<— lambdaHatCur [ i ] 
bestErro r<— error 

} 

After the loop and in case of the "flat" -mode we calculate a moving average of order two of the 
differences in the discrete curve U ^ uu and select the minimum of these values with respect to 
a penalty for big U . 

if (mcxie=" flat"){ 

diff<-diff (sigma2HatCur [1 : length ( sigma2HatCur ) ] ) 
ma<— rep (0 , length ( diff ) ) 

maOrdcr<— 2 

for(i in ( 1 : length ( diff ))) { 

ma [ i ]<— sum( abs ( d if f [max( i— maOrder , 1 ) : min( i+maOrder — 1 , length (diff 
))]))/ (min( i+maOrder , length (diff) )— max( i— maOrder , 1 ) ) 

} 

Us<-(v[2] - V [1] ) *(10+(1:( cutOfflterations -1)) ) 
alpha<— le— 5 

i<— which . min(ma+alpha*Us ) 
U<-Us [ i ] 

sigmaHat<— sqrt (pmax(0 , sigma2HatCur [ i ] ) ) 
gammaHat<— gammaHatCur [ i ] 
lambdaHat<— lambdaHatCur [ i ] 

} 

To estimate as in Q we first have to build -01^: 

z<-FT(knew-r*T, exp(-knew+r*T) *opnew ,TRUE) 
phiNu<-l-v*(v+li )*z$fy 
psiNu<— rep ( , length ( phiNu ) ) 

psiNu [( length(v)/2) : length(v) ]<-l/T*logc (phiNu [( length (v) /2) : length( 

v)]) 

psiNu [1 : ( length (v) /2-1) ]<-Re( psiNu [( length (v) -1) : ( length (v) /2+1) ] ) -1 
i*Im(psiNu [( length (v) -1) : ( length (v) /2+1) ] ) 

As before we use a loop to calculate v for different cut-off values. Since we have a separate cut-off 
UNu in the cases of "oracle" and "flat" , the number of iterations is restricted to the selected cut- 
off U of the three parameters (in that way we save calculations, because v needs smaller cut-offs 
in general). The resulting vu are stored in a matrix nuHatCur. Furthermore, we initiate some 
auxiliary variables and the f latTopKernel, as defined in (|9|. 

if (mode=" oracle" | mode="flat") 
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cut Off Iter at ions<-max( ceiling (U/ ( v [2] - v [ 1 ] ) ) ,20) -10 
nuHatCur<— matrix (O,cut0fflterations , length ( v ) ) 
nuError2< — 1 
i f (mcxie=" oracle" ) 

nuError2Cur<— rep (0 , cutOfflterations) 
if (mcxle:="PLS" ) 

objective< — 1 

flatTopKcrncl<-function(x) {(abs(x) < = 0.05) + (abs (x)<l & abs(x) >0.05)*( 

exp(-exp(-(abs(x) -0.05) ' ( -2) ) / ( abs (x) -1) ' 2) ) } 
for ( i in 1 : cutOfflterations ) { 
if (mcxie=" fix" ) 

UCur<-UfixNu 
else 

UCur <- (v[2] -V [1] ) *( i+10) 

} 

The first step in this loop is to calculate i>u with an inverse Fourier transform as in Q . Depending 
on the mode we use the final estimators a, 7 and A or the parameters per cut-off (since we have in 
the "PLS" mode a common cut-off). xCur denotes the grid on which the jump density is provided. 
It will be the same in every iteration step because it only depends on the grid v of psiNu. 

if (mode=" oracle" | mode="flat" | mode===" f ix " ) { 
lam<— lambdaHat 

fnuHat<— ( psiNu— 1 i *gammaHat*v+lam+sigmaIIat "2*v"{2}/2)* 
flatTopKernel (v/UCur) 

} 

else if (mode="PLS" ){ 
lam<— lambdaHatCur [ i ] 

fnuHat<—( psiNu— 1 i *gammaHatCur [ i ] *v+lam+sigma2IIatCur [i]*v"{2}/2)* 
flatTopKernel (v/UCur) 

} 

fxy_nu <- FT(v,fnuHat, noInverse==FALSE) 

xCur <— fxy _nu$u 

nuHatCur[i ,] <- Re( fxy _nu$ fy ) 

Up to now the condition > is not ensured and needs a correction of %. Referring to 



Belomestny and ReiQ (2006b), we seek for a ^ such that the integral /jj max{0, j>(7(a;) — £,}dx 
equals A or Xjj, respectively. We approximate the difference between the integral and lam as a 
function eq of xi, approximate its root, using the secants method, and redefine i>u as the corrected 



version. 



wNu<-rep(xCur[2] -xCur [1] , length (xCur) ) 
wNu[l]<-0.5*wNu[l] 

wNu[ length (wNu) ]<-0 . 5 *wNu [ length (wNu) ] 
if (lam>0){ 

eq<— function ( xi ) {smn(wNu* (pmax(0 , nuHatCur [ i ,] — xi ) ) )~lam} 
xiO<-0 

fxiO<-eq(xiO) 
xil<-l 

f xi 1<— cq ( xi 1 ) 
v^rhile(abs( fxil ) >0.001){ 

xi2<-xil -(xil-xiO ) / ( fxi 1 -f xi ) * fxi 1 

xiO<-xil 

xi 1<— max( , xi2 ) 
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fxiO<-fxil 
f xi 1<— eq ( xi 1 ) 

} 

} else 

xi2<-0 

nuHatCur [ i , ] <— pmax( , nuHatCur [ i ,] — xi2 ) 

Now we are in position to do the mode depending choice of i>. In the "PLS" mode we also select 
the other three parameters. In "oracle" we calculate for every UCur the L^-distance to the true v 
and pick the cut-ofF UNu with the smallest one. In "fix" -mode there is only one estimation. 

i f ( mode=" oracle" ){ 
nuTrue<— Nu(xCur ) 

nuError2Cur [ i ]<— smn(wNu* ( nuHatCur [ i ,] — nuTrue ) " 2) 
i f ( nuError2 <0 | nuError2Cur [ i ] < nuError2 ) { 

UNu<-UCur 

x<— xCur 

nuHat<— nuHatCur [ i ,] 
nuError2<— nuError2Cur [ i ] 

} 

} 

i f (mode^" fix" ) { 
UNu<-UCur 
x<— xCur 

nuHat<— nuHatCur [ i , ] 
if (model!=" real. data" ){ 
nuTrue<— Nu ( xCur ) 

nuError2<— sum(wNu* ( nuHatCur [ i ,]— nuTrue) "2) 

} 

} 

Still within the loop we deal with the "PLS" mode. Since it is driven by a least squares distance 
between the observations and the estimated model, we have to compute O from au,^u, and vu- 
It is important that the martingale condition ([T]) holds for each U , otherwise the right-hand site 
of equation ([s]) can have a singularity. Therefore, we correct Xu and vu simultaneously. Remark 
that this is not done in the other modes (and actually cannot be done before the calculation of v 



is finished). The option pricing of (kE, opE) itself is described in detail in section B.3 From all 



available strikes kE we chose the ones next to the given observations sk and get (kHat, opHat). 
if (mode=:"PLS" ){ 

beta<— (0.5 *sigma2HatCur [ i ]+gammaHatCur [ i ] ) / (lambdaHatCur [ i ] — smn( 

wNu*exp ( xCur ) *nuHatCur [ i ,]) ) 
lambdaHatCur [ i ]<— beta*lambdaHatCur [ i ] 
nuHatCur [ i , ] <— beta*nuHatCur [ i ,] 

vE< — 2'9/2 + (0:(2'9-l))*2^9/(2'9-l) 

yE!<— sapply ( vE , function (v) {ZetaHat (v , sigma2HatCur [ i ] , 

gammaHatCur [ i ] , lambdaHatCur [ i ] , xCur , nuHatCur [ i , ] ) }) 
f xy<-FT ( vE , yE , FALSE ) 
kEl<— fxy $u 
opE<-Re(fxy$fy ) 
kHat<-rep(0 , length (sk) ) 
opHat<— rep ( , length ( sk ) ) 
for(j in ( 1 : length ( sk ))) { 

index<— which . min(abs (kE—sk [ j ]) ) 

kHat [ j ]<-kE[ index] 
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opHat [ j ]<— opE [index] 

} 

} 

To penalize large second derivatives of v, the second order differences of nuHatCur are provided in 
nuDeriv. The term to be minimized is given in objectiveCur. In every iteration step we check 
whether the new objective is smaller then the last best. 

nuDeriv<-diff (nuHatCur [i ,] ,2) 

wB<— rep (xCur [2] — xCur [1] , length ( nuDeriv ) ) 

wE[l]<-0.5*wE[l] 

wE[ length (wE) ]<-0 . 5 *wE [ length (wE) ] 
alpha<-0 #le-8 

objectiveCur<— sum((opHat— snop-fpniax(0. ,1 — exp(sk— r*T) ) ) '2)+alpha* 

sum(wE* ( nuDeriv / ( xCur [2] - xCur [ 1 ] ) " 2 ) " 2 ) 
if (! is . nan( objectiveCur ) &: (objective<0 | objectiveCur <objective 

)){ 
U<-UCur 

UNu<-UCur 

sigmaHat<— sqrt (pniax(0 . sigma2HatCur [ i ] ) ) 
gammaHat<— gammaHatCur [ i ] 
lambdaHat<— lambdaHatCur [ i ] 
x<— xCur 

nuHat<— nuHatCur [ i , ] 
obj ect ivc<— obj cctiveCur 
if (model !=" real .data" ) { 
nuTrue<-Nu( xCur ) 

nuError2<— sum(wNu* ( nuHat— nuTrue ) ' 2) 

} 

} 

In the "flat" mode the choice of Ui, is made at the end of the loop. We approximate the i^-norm 
of the derivative of the map U ^ i>u and select its minimum. 

if (mode^" flat"){ 

derivativU<-diff (nuHatCur)/(v[2] -V [1]) 
derivativL2<— rep (0 , length (derivativU [ ,1]) ) 
for ( i in ( 1 : length ( derivativU [, 1 ])) ) 

derivativL2 [ i ]<—sum(wNu* derivativU [ i ,] ' 2) 
i<— which . min ( derivativL2 ) 
UNu<-(v[2]-v[l])*(i+10) 
x<— xCur 

nuHat<— nuHatCur [ i , ] 
if (modeI!=" real. data" ){ 
nuTrue<-Nu ( xCur ) 

nuError2<-sum(wNu* (nuHat— nuTrue) '2) 

} 

} 

Now we have all estimators and can return the results. 
B.3 Option pricing from Levy triplet 

Given values of cr, 7, A and v we want to provide the C-function. This happens in two situations 
of the script. On the one hand we simulate observations from a specific model and apply the 
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calibration to the noised data and on the other hand, we have to calculate option prices from the 



estimated parameters in the least squares method. In (Cont and Tankov 2004b p. 361 et seq.) a 
good description of this option pricing can be found. 
We define the function 

Civ) :=e--MLzl^ 
tv(l + IV) 

which is implemented as Zeta(v) and ZetaHat (v, sigma2Hat, gaiimiaHat, lambdaHat, x, nuHat), 

respectively. Then O as function of log-strike k is obtained by inverse Fourier transform of (. 
Hence, we evaluate C on a discrete grid x to get pairs of log-strikes and option prices (k, op). 
Because of the FFT algorithm a finer grid x leads to a larger limits of k and the other way around 
a higher range of x implies a smaller mesh of k. 

A<-2~10 
MC-12 
]N<-2"M 
K-0: (N-1) 
Delta<-A/(N-1) 
x< — A/2+l*Delta 
y<— Zeta (x) 
fxy<-FT(x,y, FALSE) 
k<— fxy $u 

op<-Re( fxySfy )4pmax(0. , 1 -exp (k-r*T) ) 

For simulations we noise and sample all pairs (k,op) and apply the calibration procedure in every 
Monte Carlo iteration step. Here the sampling is done according to a normal distribution centered 
at at-the-money strikes. Depending on design, the strikes a distributed deterministic or random. 

for (mi in 1 : monteCarlo ) { 

nop<— op+rnorm ( length (op) ,0 , ( noiseLevel*abs (Rb( fxy$fy ) ) ) ) 



Ml<— round ( length ( op ) / 2 )— which . min( abs ( op [ 1 : round ( length ( op ) / 2) 

]-10"(-6))) 
M2<-which . niax( op ) +M1 
MK-which . niax( op ) -Ml 

i f ( design=" random" ) { 

prob<-exp(-(k [M2:Ml]-r*T) "2) 
prob<— prob /sum( prob ) 

ind _sub<— sort (sample (M2:M1, sampleSize , prob=prob ) ) 
sk<— k [ ind _sub ] 
snop<— nop [ ind _sub] 
}else if(desig n=" deterministic" ){ 

quant ils<—qnorm( c (l:sampleSize)/(sampleSize+l),0,2"(— 0.5)) 
ind _sub<— sapply ( quant ils , function (x) {which . min( abs (k— r*T— x) ) }) 
sk<— k [ ind _sub ] 
snop<— nop [ ind _sub ] 

} 

snop<— snop4pniax(0 . ,1 — exp(sk— r*T) ) 

} 
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